## -------------------------------------------------------------------
## Firm-level analysis 
## -------------------------------------------------------------------
rm(list=ls())

## library
library(sandwich)
library(dplyr)
library(gridExtra)
library(ggplot2)

## setup -------------------------------------------------------------
## setting the working directory as the current directory
setwd(getwd())

## check/create the folder to put outputs
if (file.exists(file.path(getwd(), "figs"))){
    FIG_DIR <- file.path(getwd(), "figs")
} else {
    dir.create(file.path(file.path(getwd(), "figs")))
    FIG_DIR <- file.path(getwd(), "figs")
}
    
source("util.R")

resptwice <- c("R_3lKP5vVopq27Q5T", "R_ewwjPKYf1mQtVkx",
               "R_8wVpt534HFyvbG5","R_7QwtntdSZ7ONCaV",
               "R_787WpeOkIatHWkZ", "R_8umeWuf3qISiJcV",
               "R_85MeGXniJgeyA6h", "R_0Sy2CDRAw2ddHql",
               "R_5pDS87UuBf5ekU5","R_9YwTFMe01hCQIDP",
               "R_38WL2rp1ocK9c1v","R_0qEMbRZcn0W3gMd",
               "R_3UdWL2nJMFa61hP","R_ctKLc9LF3hANh4x")

## read survey data: the original survey data
survey <- read.csv(file.path("crsurvey_final.csv"),
                   stringsAsFactors=FALSE)
## check
which(survey$V1 %in% resptwice)

Data <- read.csv(file.path("conjoint_final.csv"),
                   stringsAsFactors=FALSE)


## setting levels
Data$dispute <- factor(Data$dispute, levels=c("passive", "moderate", "aggresive"))
Data$flex <- factor(Data$flex, levels=c("low", "moderate", "high"))
Data$reduc_barrier <- factor(Data$reduc_barrier, levels=c("low", "moderate", "high"))
Data$invest <- factor(Data$invest, levels=c("weak", "moderate", "strong"))
Data$subsidy <- factor(Data$subsidy, levels=c("low", "moderate", "high"))

## running separate regression
varnames <- c("invest", "reduc_barrier", "subsidy",
              "dispute", "flex")
baselines <- c("weak", "low", "low", "passive", "low")


## firms in tradable industry AND responded to the conjoint: 214 firms
resp <- unique(Data$resp_id)
ids.final <- which(resp %in% survey$V1)
length(ids.final)

ids <- resp[ids.final]


survey2 <- survey[which(survey$V1 %in%ids),]


## -------------------------------------------------------------------
## Figure 5: Domestic vs. Exporter/Multinational Difference
## -------------------------------------------------------------------

diff <- TRUE
flip <- TRUE ## exporter/multinational - domestic

title <- c("Domestic",
           "Exporter/Multinational",
           "Difference")

filename <- "figure5.pdf"


idx1 <- which(survey2$iskcat2 == "Domestic")
idx21 <- which(survey2$iskcat2=="Exporter")
idx22 <- which((survey2$iskcat2=="Exporter" & survey2$own_3 > 80) | survey2$iskcat2=="Multinational")
idx2 <- c(idx21,idx22)


## creating Figure 5
length(idx1)
length(idx2)
if(!diff){
    length(idx3)
}


resp_id1 <- as.character(survey2$V1[idx1])
resp_id2 <- as.character(survey2$V1[idx2])
if(!diff){
    resp_id3 <- as.character(survey2$V1[idx3])
}

Data1 <- Data[Data$resp_id%in%resp_id1,]
Data2 <- Data[Data$resp_id%in%resp_id2,]
if(!diff){
    Data3 <- Data[Data$resp_id%in%resp_id3,]
}

length(unique(Data1$resp_id))
length(unique(Data2$resp_id))
if(!diff){
    length(unique(Data3$resp_id))
}

if(!diff){
    counts <- c(nrow(Data1)/2, nrow(Data2)/2,
                nrow(Data3)/2)
} else {
    counts <- c(nrow(Data1)/2, nrow(Data2)/2,
                (nrow(Data1)+nrow(Data2))/2)

}

main.cex <- 1.2
title <- paste(title, "\n (", counts, ")", sep="")


results <- list()

if(!diff){
    for(j in 1:3){
        cmd <- paste("df <- Data", j, sep="")
        eval(parse(text= cmd))
        cmd <- paste("df$resp_id <- as.factor(as.character(df$resp_id))", sep="")
        eval(parse(text= cmd))

        for(i in 1:length(varnames)){
            formula <- paste("select ~ ", varnames[i], sep="")
            mod <- lm(as.formula(formula), data=df)
            pe.i <- as.matrix(mod$coef[-1])
            base.i <- matrix(0, nrow=1, ncol=2)
            rownames(base.i) <- paste(varnames[i], baselines[i], sep="")

                                        # check: observations are too few for cluster se.
            if(length(unique(df$resp_id))>50){
                vcov_mat <- cluster_se(mod, df$resp_id)
                se.i <- as.matrix(sqrt(diag(vcov_mat))[-1])
            } else {
                se.i <- as.matrix(sqrt(diag(vcov(mod)))[-1])
            }

            ## combining results
            result.i <- cbind(pe.i, se.i)
            result.i <- rbind(base.i, result.i)
            if(i==1){
                result <- result.i
            } else {
                result <- rbind(result, result.i)
            }
        }
        results[[j]] <- result
    }
} else {
    for(j in 1:2){
        cmd <- paste("df <- Data", j, sep="")
        eval(parse(text= cmd))
        cmd <- paste("df$resp_id <- as.factor(as.character(df$resp_id))", sep="")
        eval(parse(text= cmd))

        for(i in 1:length(varnames)){
            formula <- paste("select ~ ", varnames[i], sep="")
            mod <- lm(as.formula(formula), data=df)
            pe.i <- as.matrix(mod$coef[-1])
            base.i <- matrix(0, nrow=1, ncol=2)
            rownames(base.i) <- paste(varnames[i], baselines[i], sep="")

            ## check: observations are too few for cluster se.
            if(length(unique(df$resp_id))>50){
                vcov_mat <- cluster_se(mod, df$resp_id)
                se.i <- as.matrix(sqrt(diag(vcov_mat))[-1])
            } else {
                se.i <- as.matrix(sqrt(diag(vcov(mod)))[-1])
            }

            ## combining results
            result.i <- cbind(pe.i, se.i)
            result.i <- rbind(base.i, result.i)
            if(i==1){
                result <- result.i
            } else {
                result <- rbind(result, result.i)
            }
        }
        results[[j]] <- result
    }

    
    Data_diff <- rbind(Data1, Data2)
    Data_diff$resp_id <- as.factor(as.character(Data_diff$resp_id))
    Data_diff$varcol <- c(rep(1, nrow(Data1)), rep(0, nrow(Data2)))
    resp.idx <- which(colnames(survey) == "V1")

    ## regression diff
    for(i in 1:length(varnames)){
        formula <- paste("select ~ ", varnames[i], "+ varcol +", varnames[i], ":varcol", sep="")
        mod <- lm(as.formula(formula), data=Data_diff)
        pe.i <- as.matrix(mod$coef[-(1:4)])
        ## flip the order
        if(flip){
            pe.i <- -1*pe.i
        }
        base.i <- matrix(0, nrow=1, ncol=2)
        rownames(base.i) <- paste(varnames[i], baselines[i], sep="")
        vcov_mat <- cluster_se(mod, Data_diff$resp_id)[-(1:4),-(1:4)]
        se.i <- as.matrix(sqrt(diag(vcov_mat)))
        result.i <- cbind(pe.i, se.i)
        result.i <- rbind(base.i, result.i)
        if(i==1){
            result3 <- result.i
        } else {
            result3 <- rbind(result3, result.i)
        }
        
    }
    results[[3]] <- result3
}



## creating graphs ------------------------------------------

main.cex=2.5
col1 <- "black"
pdf(file.path(FIG_DIR, filename),
    width=22,height=10) # for paper

par(mfrow=c(1,4))

for(p in 0:length(results)){

    if(p==0){
        par(cex.lab=1.3,
            mar=c(5,28,6,2.5),
            cex.lab=1.7, cex.axis=1.7,
            cex.main=main.cex)
        
        yloc <- c(23,22,21,18,17,16,13,12,11,8,7,6,3,2,1)
        pe <-  as.numeric(results[[1]][,1])
        plot(pe,yloc,
             ylim=c(1,24),
             xlim=c(0,0.001),
             xaxt='n',yaxt='n',bty='n',pch='',ylab='',xlab='')

        
        ## investment Provision
        mtext("Investment Protection", side=2, font=2,
              las=1, line=25, cex=1.75, at=24, adj=0, col=col1)
        mtext("Weak Legal Protection", side=2,
              las=1, line=15, cex=1.5, at=23, adj=0)
        mtext("Moderate Legal Protection", side=2,
              las=1, line=15, cex=1.5, at=22, adj=0)
        mtext("Strong Legal Protection", side=2,
              las=1, line=15, cex=1.5, at=21, adj=0)

        ## Reduction of tariff Barriers
        mtext("Reduction of Trade Barriers", side=2, font=2,
              las=1, line=25, cex=1.75, at=19, adj=0, col=col1)
        mtext("Low Reduction", side=2,
              las=1, line=15, cex=1.5, at=18, adj=0)
        mtext("Moderate Reduction", side=2,
              las=1, line=15, cex=1.5, at=17, adj=0)
        mtext("High Reduction", side=2,
              las=1, line=15, cex=1.5, at=16, adj=0)

        ## Export Subsidies
        mtext("Export Subsidies", side=2, font=2,
              las=1, line=25, cex=1.75, at=14, adj=0, col=col1)
        mtext("Low Subsidies", side=2,
              las=1, line=15, cex=1.5, at=13, adj=0)
        mtext("Moderate Subsidies", side=2,
              las=1, line=15, cex=1.5, at=12, adj=0)
        mtext("High Subsidies", side=2,
              las=1, line=15, cex=1.5, at=11, adj=0)

        ## Dispute Settlement
        mtext("Dispute Settlement Mechanism", side=2, font=2,
              las=1, line=25, cex=1.75, at=9, adj=0, col=col1)
        mtext("Passive Use", side=2,
              las=1, line=15, cex=1.5, at=8, adj=0)
        mtext("Moderate Use", side=2,
              las=1, line=15, cex=1.5, at=7, adj=0)
        mtext("Aggressive Use", side=2,
              las=1, line=15, cex=1.5, at=6, adj=0)

        ## Dispute Settlement
        mtext("Flexibility of Commitments", side=2, font=2,
              las=1, line=25, cex=1.75, at=4, adj=0, col=col1)
        mtext("Low Flexibility", side=2,
              las=1, line=15, cex=1.5, at=3, adj=0)
        mtext("Moderate Flxibiilty", side=2,
              las=1, line=15, cex=1.5, at=2, adj=0)
        mtext("High Flexibility", side=2,
              las=1, line=15, cex=1.5, at=1, adj=0)

    } else {
        par(cex.lab=1.3,
            mar=c(6,4,6,2.5),
            cex.lab=2, cex.axis=2,
            cex.main=main.cex)

        ## yloc <- seq(15, by=-1, length.out=15)
        yloc <- c(23,22,21,18,17,16,13,12,11,8,7,6,3,2,1)
        col1= "black"
        pe <-  as.numeric(results[[p]][,1])
        
        lo <- as.numeric(results[[p]][,1] - qnorm(.975)*results[[p]][,2])
        hi <- as.numeric(results[[p]][,1] + qnorm(.975)*results[[p]][,2])
        
        colors <- rep("black", length(pe))

        colors[c(3,12)] <- "red"
        ## colors[c(6)] <- "red"
        
        plot(pe, yloc,
             xlim=c(-0.4, 0.5),
             ylim=c(1,24),
             pch=16,
             main=title[p],
             yaxt="n", ylab="", cex=2.5,
             col=colors,
             xlab="",
             )

        mtext("Change in Pr(Support Policy)",
              side=1, 
              line=4.5, cex=1.65)
        
        count <- 1
        for(i in 1:length(yloc)){
            loc <- yloc[i]
            segments(lo[count], loc, hi[count], loc, lwd=2, col=colors[i])
            count <- count + 1
        }
        abline(v=0, lty=2)
        xloc <- seq(-0.4,0.5,by=.1)
        for(k in xloc){
            abline(v=k, lty=2, col="grey80")
        }

        axis(side=2, at=yloc, tick=TRUE, labels=rep("",length(yloc)))
             

    }
}
dev.off()



## -------------------------------------------------------------------
## Figure 7: Net exporting vs Net importing
## -------------------------------------------------------------------

diff <- TRUE
flip <- FALSE # exporting - importing

title <- c("Net Exporting Industry",
           "Net Importing Industry",
           "Difference")

filename <- "figure7.pdf"

net <- survey2$main_exp - survey2$main_imp
total <- survey2$main_exp + survey2$main_imp
idx1 <- which(total > 50000000 & net > 0)
idx2 <- which(total > 50000000 & net < 0)


## creating Figure 7
length(idx1)
length(idx2)
if(!diff){
    length(idx3)
}


resp_id1 <- as.character(survey2$V1[idx1])
resp_id2 <- as.character(survey2$V1[idx2])
if(!diff){
    resp_id3 <- as.character(survey2$V1[idx3])
}

Data1 <- Data[Data$resp_id%in%resp_id1,]
Data2 <- Data[Data$resp_id%in%resp_id2,]
if(!diff){
    Data3 <- Data[Data$resp_id%in%resp_id3,]
}

length(unique(Data1$resp_id))
length(unique(Data2$resp_id))
if(!diff){
    length(unique(Data3$resp_id))
}

if(!diff){
    counts <- c(nrow(Data1)/2, nrow(Data2)/2,
                nrow(Data3)/2)
} else {
    counts <- c(nrow(Data1)/2, nrow(Data2)/2,
                (nrow(Data1)+nrow(Data2))/2)

}

main.cex <- 1.2
title <- paste(title, "\n (", counts, ")", sep="")


results <- list()

if(!diff){
    for(j in 1:3){
        cmd <- paste("df <- Data", j, sep="")
        eval(parse(text= cmd))
        cmd <- paste("df$resp_id <- as.factor(as.character(df$resp_id))", sep="")
        eval(parse(text= cmd))

        for(i in 1:length(varnames)){
            formula <- paste("select ~ ", varnames[i], sep="")
            mod <- lm(as.formula(formula), data=df)
            pe.i <- as.matrix(mod$coef[-1])
            base.i <- matrix(0, nrow=1, ncol=2)
            rownames(base.i) <- paste(varnames[i], baselines[i], sep="")

            ## check: observations are too few for cluster se.
            if(length(unique(df$resp_id))>50){
                vcov_mat <- cluster_se(mod, df$resp_id)
                se.i <- as.matrix(sqrt(diag(vcov_mat))[-1])
            } else {
                se.i <- as.matrix(sqrt(diag(vcov(mod)))[-1])
            }

            ## combining results
            result.i <- cbind(pe.i, se.i)
            result.i <- rbind(base.i, result.i)
            if(i==1){
                result <- result.i
            } else {
                result <- rbind(result, result.i)
            }
        }
        results[[j]] <- result
    }
} else {
    for(j in 1:2){
        cmd <- paste("df <- Data", j, sep="")
        eval(parse(text= cmd))
        cmd <- paste("df$resp_id <- as.factor(as.character(df$resp_id))", sep="")
        eval(parse(text= cmd))

        for(i in 1:length(varnames)){
            formula <- paste("select ~ ", varnames[i], sep="")
            mod <- lm(as.formula(formula), data=df)
            pe.i <- as.matrix(mod$coef[-1])
            base.i <- matrix(0, nrow=1, ncol=2)
            rownames(base.i) <- paste(varnames[i], baselines[i], sep="")

                                        # check: observations are too few for cluster se.
            if(length(unique(df$resp_id))>50){
                vcov_mat <- cluster_se(mod, df$resp_id)
                se.i <- as.matrix(sqrt(diag(vcov_mat))[-1])
            } else {
                se.i <- as.matrix(sqrt(diag(vcov(mod)))[-1])
            }

            ## combining results
            result.i <- cbind(pe.i, se.i)
            result.i <- rbind(base.i, result.i)
            if(i==1){
                result <- result.i
            } else {
                result <- rbind(result, result.i)
            }
        }
        results[[j]] <- result
    }

    
    Data_diff <- rbind(Data1, Data2)
    Data_diff$resp_id <- as.factor(as.character(Data_diff$resp_id))
    Data_diff$varcol <- c(rep(1, nrow(Data1)), rep(0, nrow(Data2)))
    resp.idx <- which(colnames(survey) == "V1")

    ## regression diff
    for(i in 1:length(varnames)){
        formula <- paste("select ~ ", varnames[i], "+ varcol +", varnames[i], ":varcol", sep="")
        mod <- lm(as.formula(formula), data=Data_diff)
        pe.i <- as.matrix(mod$coef[-(1:4)])
        ## flip the order
        if(flip){
            pe.i <- -1*pe.i
        }
        base.i <- matrix(0, nrow=1, ncol=2)
        rownames(base.i) <- paste(varnames[i], baselines[i], sep="")
        vcov_mat <- cluster_se(mod, Data_diff$resp_id)[-(1:4),-(1:4)]
        se.i <- as.matrix(sqrt(diag(vcov_mat)))
        result.i <- cbind(pe.i, se.i)
        result.i <- rbind(base.i, result.i)
        if(i==1){
            result3 <- result.i
        } else {
            result3 <- rbind(result3, result.i)
        }
        
    }
    results[[3]] <- result3
}



## creating graphs ------------------------------------------

main.cex=2.5
col1 <- "black"
pdf(file.path(FIG_DIR, filename),
    width=22,height=10) # for paper

par(mfrow=c(1,4))

for(p in 0:length(results)){

    if(p==0){
        par(cex.lab=1.3,
            mar=c(5,28,6,2.5),
            cex.lab=1.7, cex.axis=1.7,
            cex.main=main.cex)
        
        yloc <- c(23,22,21,18,17,16,13,12,11,8,7,6,3,2,1)
        pe <-  as.numeric(results[[1]][,1])
        plot(pe,yloc,
             ylim=c(1,24),
             xlim=c(0,0.001),
             xaxt='n',yaxt='n',bty='n',pch='',ylab='',xlab='')

        
        ## investment Provision
        mtext("Investment Protection", side=2, font=2,
              las=1, line=25, cex=1.75, at=24, adj=0, col=col1)
        mtext("Weak Legal Protection", side=2,
              las=1, line=15, cex=1.5, at=23, adj=0)
        mtext("Moderate Legal Protection", side=2,
              las=1, line=15, cex=1.5, at=22, adj=0)
        mtext("Strong Legal Protection", side=2,
              las=1, line=15, cex=1.5, at=21, adj=0)

        ## Reduction of tariff Barriers
        mtext("Reduction of Trade Barriers", side=2, font=2,
              las=1, line=25, cex=1.75, at=19, adj=0, col=col1)
        mtext("Low Reduction", side=2,
              las=1, line=15, cex=1.5, at=18, adj=0)
        mtext("Moderate Reduction", side=2,
              las=1, line=15, cex=1.5, at=17, adj=0)
        mtext("High Reduction", side=2,
              las=1, line=15, cex=1.5, at=16, adj=0)

        ## Export Subsidies
        mtext("Export Subsidies", side=2, font=2,
              las=1, line=25, cex=1.75, at=14, adj=0, col=col1)
        mtext("Low Subsidies", side=2,
              las=1, line=15, cex=1.5, at=13, adj=0)
        mtext("Moderate Subsidies", side=2,
              las=1, line=15, cex=1.5, at=12, adj=0)
        mtext("High Subsidies", side=2,
              las=1, line=15, cex=1.5, at=11, adj=0)

        ## Dispute Settlement
        mtext("Dispute Settlement Mechanism", side=2, font=2,
              las=1, line=25, cex=1.75, at=9, adj=0, col=col1)
        mtext("Passive Use", side=2,
              las=1, line=15, cex=1.5, at=8, adj=0)
        mtext("Moderate Use", side=2,
              las=1, line=15, cex=1.5, at=7, adj=0)
        mtext("Aggressive Use", side=2,
              las=1, line=15, cex=1.5, at=6, adj=0)

        ## Dispute Settlement
        mtext("Flexibility of Commitments", side=2, font=2,
              las=1, line=25, cex=1.75, at=4, adj=0, col=col1)
        mtext("Low Flexibility", side=2,
              las=1, line=15, cex=1.5, at=3, adj=0)
        mtext("Moderate Flxibiilty", side=2,
              las=1, line=15, cex=1.5, at=2, adj=0)
        mtext("High Flexibility", side=2,
              las=1, line=15, cex=1.5, at=1, adj=0)

    } else {
        par(cex.lab=1.3,
            mar=c(6,4,6,2.5),
            cex.lab=2, cex.axis=2,
            cex.main=main.cex)

        ## yloc <- seq(15, by=-1, length.out=15)
        yloc <- c(23,22,21,18,17,16,13,12,11,8,7,6,3,2,1)
        col1= "black"
        pe <-  as.numeric(results[[p]][,1])
        
        lo <- as.numeric(results[[p]][,1] - qnorm(.975)*results[[p]][,2])
        hi <- as.numeric(results[[p]][,1] + qnorm(.975)*results[[p]][,2])
        
        colors <- rep("black", length(pe))

        plot(pe, yloc,
             xlim=c(-0.4, 0.5),
             ylim=c(1,24),
             pch=16,
             main=title[p],
             yaxt="n", ylab="", cex=2.5,
             col=colors,
             xlab="",
             )

        mtext("Change in Pr(Support Policy)",
              side=1, 
              line=4.5, cex=1.65)
        
        count <- 1
        for(i in 1:length(yloc)){
            loc <- yloc[i]
            segments(lo[count], loc, hi[count], loc, lwd=2, col=colors[i])
            count <- count + 1
        }
        abline(v=0, lty=2)
        xloc <- seq(-0.4,0.5,by=.1)
        for(k in xloc){
            abline(v=k, lty=2, col="grey80")
        }

        axis(side=2, at=yloc, tick=TRUE, labels=rep("",length(yloc)))
             

    }
}
dev.off()



## -------------------------------------------------------------------
## Figure 8: Agriculture vs Manufacturing industries
## -------------------------------------------------------------------


procomer <<- read.csv(file.path("procomerALL.csv"))
isic2 <- mapply(getISIC2, survey2$isic)
tradables <- c(paste("0", 1:9, sep=""), as.character(10:33))


diff <- TRUE

title <- c("Agriculture Industry",
           "Manufacturing Industry",
           "Difference")

filename <- "figure8.pdf"

agrExp <- as.numeric(mapply(IsAgrExporter, survey2$proccode))
agric <- c(paste("0", 1:3, sep=""))
manuf <- as.character(10:33)
idx1 <- which((is.na(agrExp) & isic2 %in% agric) | agrExp==1)
idx2 <- which((is.na(agrExp) & isic2 %in% manuf) | agrExp==0)



## creating Figure 8
length(idx1)
length(idx2)
if(!diff){
    length(idx3)
}


resp_id1 <- as.character(survey2$V1[idx1])
resp_id2 <- as.character(survey2$V1[idx2])
if(!diff){
    resp_id3 <- as.character(survey2$V1[idx3])
}

Data1 <- Data[Data$resp_id%in%resp_id1,]
Data2 <- Data[Data$resp_id%in%resp_id2,]
if(!diff){
    Data3 <- Data[Data$resp_id%in%resp_id3,]
}

length(unique(Data1$resp_id))
length(unique(Data2$resp_id))
if(!diff){
    length(unique(Data3$resp_id))
}

if(!diff){
    counts <- c(nrow(Data1)/2, nrow(Data2)/2,
                nrow(Data3)/2)
} else {
    counts <- c(nrow(Data1)/2, nrow(Data2)/2,
                (nrow(Data1)+nrow(Data2))/2)

}

main.cex <- 1.2
title <- paste(title, "\n (", counts, ")", sep="")


results <- list()

if(!diff){
    for(j in 1:3){
        cmd <- paste("df <- Data", j, sep="")
        eval(parse(text= cmd))
        cmd <- paste("df$resp_id <- as.factor(as.character(df$resp_id))", sep="")
        eval(parse(text= cmd))

        for(i in 1:length(varnames)){
            formula <- paste("select ~ ", varnames[i], sep="")
            mod <- lm(as.formula(formula), data=df)
            pe.i <- as.matrix(mod$coef[-1])
            base.i <- matrix(0, nrow=1, ncol=2)
            rownames(base.i) <- paste(varnames[i], baselines[i], sep="")

            ## check: observations are too few for cluster se.
            if(length(unique(df$resp_id))>50){
                vcov_mat <- cluster_se(mod, df$resp_id)
                se.i <- as.matrix(sqrt(diag(vcov_mat))[-1])
            } else {
                se.i <- as.matrix(sqrt(diag(vcov(mod)))[-1])
            }

            ## combining results
            result.i <- cbind(pe.i, se.i)
            result.i <- rbind(base.i, result.i)
            if(i==1){
                result <- result.i
            } else {
                result <- rbind(result, result.i)
            }
        }
        results[[j]] <- result
    }
} else {
    for(j in 1:2){
        cmd <- paste("df <- Data", j, sep="")
        eval(parse(text= cmd))
        cmd <- paste("df$resp_id <- as.factor(as.character(df$resp_id))", sep="")
        eval(parse(text= cmd))

        for(i in 1:length(varnames)){
            formula <- paste("select ~ ", varnames[i], sep="")
            mod <- lm(as.formula(formula), data=df)
            pe.i <- as.matrix(mod$coef[-1])
            base.i <- matrix(0, nrow=1, ncol=2)
            rownames(base.i) <- paste(varnames[i], baselines[i], sep="")

            ## check: observations are too few for cluster se.
            if(length(unique(df$resp_id))>50){
                vcov_mat <- cluster_se(mod, df$resp_id)
                se.i <- as.matrix(sqrt(diag(vcov_mat))[-1])
            } else {
                se.i <- as.matrix(sqrt(diag(vcov(mod)))[-1])
            }

            ## combining results
            result.i <- cbind(pe.i, se.i)
            result.i <- rbind(base.i, result.i)
            if(i==1){
                result <- result.i
            } else {
                result <- rbind(result, result.i)
            }
        }
        results[[j]] <- result
    }

    
    Data_diff <- rbind(Data1, Data2)
    Data_diff$resp_id <- as.factor(as.character(Data_diff$resp_id))
    Data_diff$varcol <- c(rep(1, nrow(Data1)), rep(0, nrow(Data2)))
    resp.idx <- which(colnames(survey) == "V1")

    ## regression diff
    for(i in 1:length(varnames)){
        formula <- paste("select ~ ", varnames[i], "+ varcol +", varnames[i], ":varcol", sep="")
        mod <- lm(as.formula(formula), data=Data_diff)
        pe.i <- as.matrix(mod$coef[-(1:4)])
        ## flip the order
        if(flip){
            pe.i <- -1*pe.i
        }
        base.i <- matrix(0, nrow=1, ncol=2)
        rownames(base.i) <- paste(varnames[i], baselines[i], sep="")
        vcov_mat <- cluster_se(mod, Data_diff$resp_id)[-(1:4),-(1:4)]
        se.i <- as.matrix(sqrt(diag(vcov_mat)))
        result.i <- cbind(pe.i, se.i)
        result.i <- rbind(base.i, result.i)
        if(i==1){
            result3 <- result.i
        } else {
            result3 <- rbind(result3, result.i)
        }
        
    }
    results[[3]] <- result3
}



## creating graphs ------------------------------------------

main.cex=2.5
col1 <- "black"
pdf(file.path(FIG_DIR, filename),
    width=22,height=10) # for paper

par(mfrow=c(1,4))

for(p in 0:length(results)){

    if(p==0){
        par(cex.lab=1.3,
            mar=c(5,28,6,2.5),
            cex.lab=1.7, cex.axis=1.7,
            cex.main=main.cex)
        
        yloc <- c(23,22,21,18,17,16,13,12,11,8,7,6,3,2,1)
        pe <-  as.numeric(results[[1]][,1])
        plot(pe,yloc,
             ylim=c(1,24),
             xlim=c(0,0.001),
             xaxt='n',yaxt='n',bty='n',pch='',ylab='',xlab='')

        
        ## investment Provision
        mtext("Investment Protection", side=2, font=2,
              las=1, line=25, cex=1.75, at=24, adj=0, col=col1)
        mtext("Weak Legal Protection", side=2,
              las=1, line=15, cex=1.5, at=23, adj=0)
        mtext("Moderate Legal Protection", side=2,
              las=1, line=15, cex=1.5, at=22, adj=0)
        mtext("Strong Legal Protection", side=2,
              las=1, line=15, cex=1.5, at=21, adj=0)

        ## Reduction of tariff Barriers
        mtext("Reduction of Trade Barriers", side=2, font=2,
              las=1, line=25, cex=1.75, at=19, adj=0, col=col1)
        mtext("Low Reduction", side=2,
              las=1, line=15, cex=1.5, at=18, adj=0)
        mtext("Moderate Reduction", side=2,
              las=1, line=15, cex=1.5, at=17, adj=0)
        mtext("High Reduction", side=2,
              las=1, line=15, cex=1.5, at=16, adj=0)

        ## Export Subsidies
        mtext("Export Subsidies", side=2, font=2,
              las=1, line=25, cex=1.75, at=14, adj=0, col=col1)
        mtext("Low Subsidies", side=2,
              las=1, line=15, cex=1.5, at=13, adj=0)
        mtext("Moderate Subsidies", side=2,
              las=1, line=15, cex=1.5, at=12, adj=0)
        mtext("High Subsidies", side=2,
              las=1, line=15, cex=1.5, at=11, adj=0)

        ## Dispute Settlement
        mtext("Dispute Settlement Mechanism", side=2, font=2,
              las=1, line=25, cex=1.75, at=9, adj=0, col=col1)
        mtext("Passive Use", side=2,
              las=1, line=15, cex=1.5, at=8, adj=0)
        mtext("Moderate Use", side=2,
              las=1, line=15, cex=1.5, at=7, adj=0)
        mtext("Aggressive Use", side=2,
              las=1, line=15, cex=1.5, at=6, adj=0)

        ## Dispute Settlement
        mtext("Flexibility of Commitments", side=2, font=2,
              las=1, line=25, cex=1.75, at=4, adj=0, col=col1)
        mtext("Low Flexibility", side=2,
              las=1, line=15, cex=1.5, at=3, adj=0)
        mtext("Moderate Flxibiilty", side=2,
              las=1, line=15, cex=1.5, at=2, adj=0)
        mtext("High Flexibility", side=2,
              las=1, line=15, cex=1.5, at=1, adj=0)

    } else {
        par(cex.lab=1.3,
            mar=c(6,4,6,2.5),
            cex.lab=2, cex.axis=2,
            cex.main=main.cex)

        ## yloc <- seq(15, by=-1, length.out=15)
        yloc <- c(23,22,21,18,17,16,13,12,11,8,7,6,3,2,1)
        col1= "black"
        pe <-  as.numeric(results[[p]][,1])
        
        lo <- as.numeric(results[[p]][,1] - qnorm(.975)*results[[p]][,2])
        hi <- as.numeric(results[[p]][,1] + qnorm(.975)*results[[p]][,2])
        
        colors <- rep("black", length(pe))

        plot(pe, yloc,
             xlim=c(-0.4, 0.5),
             ylim=c(1,24),
             pch=16,
             main=title[p],
             yaxt="n", ylab="", cex=2.5,
             col=colors,
             xlab="",
             )

        mtext("Change in Pr(Support Policy)",
              side=1, 
              line=4.5, cex=1.65)
        
        count <- 1
        for(i in 1:length(yloc)){
            loc <- yloc[i]
            segments(lo[count], loc, hi[count], loc, lwd=2, col=colors[i])
            count <- count + 1
        }
        abline(v=0, lty=2)
        xloc <- seq(-0.4,0.5,by=.1)
        for(k in xloc){
            abline(v=k, lty=2, col="grey80")
        }

        axis(side=2, at=yloc, tick=TRUE, labels=rep("",length(yloc)))
             

    }
}
dev.off()


## -------------------------------------------------------------------
## Figure 6
## -------------------------------------------------------------------

infoD <- read.csv(file.path("procomer_resp_info.csv"))

survey$V1[which(survey$V1 == "E2616; E4532")] <- "E2616"
survey1 <- merge(survey, infoD, by.x="proccode", by.y="procomerID",
                 all.x=T)

## firms in tradable industry AND responded to the conjoint: 214 firms
resp <- unique(Data$resp_id)
ids.final <- which(resp %in% survey$V1)
length(ids.final)
ids <- resp[ids.final]

survey2 <- survey[which(survey$V1 %in%ids),]


## Exporter import and export simultaneously ----------------

title <- c("Domestic",
           "Autonomous Exporter",           
           "Exporter in Global Value Chains",
           "Multinational")

filename <- "figure6.pdf"

idx1 <- which(survey2$iskcat2 == "Domestic")
idx2 <- which(survey2$iskcat2=="Exporter" & survey2$import ==2)
idx3 <- which(survey2$iskcat2=="Exporter" & survey2$import ==1)
idx4 <- which((survey2$iskcat2=="Exporter" & survey2$own_3 > 80) | survey2$iskcat2=="Multinational")


## make sure that there is no overlap
id2.id3.overlap <- which(idx2 %in% idx3)
if(length(id2.id3.overlap)>0){
    idx2 <- idx2[-id2.id3.overlap]
}

id1.id4.overlap <- which(idx1 %in% idx4)
if(length(id1.id4.overlap)>0){
    idx1 <- idx1[-id1.id4.overlap]
}

id2.id4.overlap <- which(idx2 %in% idx4)
if(length(id2.id4.overlap)>0){
    idx2 <- idx2[-id2.id4.overlap]
}

id3.id4.overlap <- which(idx3 %in% idx4)
if(length(id3.id4.overlap)>0){
    idx3 <- idx3[-id3.id4.overlap]
}


length(idx1)
length(idx2)
length(idx3)
length(idx4)

resp_id1 <- as.character(survey2$V1[idx1])
resp_id2 <- as.character(survey2$V1[idx2])
resp_id3 <- as.character(survey2$V1[idx3])
resp_id4 <- as.character(survey2$V1[idx4])

Data1 <- Data[Data$resp_id%in%resp_id1,]
Data2 <- Data[Data$resp_id%in%resp_id2,]
Data3 <- Data[Data$resp_id%in%resp_id3,]
Data4 <- Data[Data$resp_id%in%resp_id4,]

length(unique(Data1$resp_id))
length(unique(Data2$resp_id))
length(unique(Data3$resp_id))
length(unique(Data4$resp_id))


counts <- c(nrow(Data1)/2, nrow(Data2)/2,
            nrow(Data3)/2, nrow(Data4)/2)            

main.cex <- 1.2
title <- paste(title, "\n (", counts, ")", sep="")


results <- list()

for(j in 1:4){
    cmd <- paste("df <- Data", j, sep="")
    eval(parse(text= cmd))
    cmd <- paste("df$resp_id <- as.factor(as.character(df$resp_id))", sep="")
    eval(parse(text= cmd))

    for(i in 1:length(varnames)){
        formula <- paste("select ~ ", varnames[i], sep="")
        mod <- lm(as.formula(formula), data=df)
        pe.i <- as.matrix(mod$coef[-1])
        base.i <- matrix(0, nrow=1, ncol=2)
        rownames(base.i) <- paste(varnames[i], baselines[i], sep="")

        # check: observations are too few for cluster se.
        if(length(unique(df$resp_id))>50){
            vcov_mat <- cluster_se(mod, df$resp_id)
            se.i <- as.matrix(sqrt(diag(vcov_mat))[-1])
        } else {
            se.i <- as.matrix(sqrt(diag(vcov(mod)))[-1])
        }

        ## combining results
        result.i <- cbind(pe.i, se.i)
        result.i <- rbind(base.i, result.i)
        if(i==1){
            result <- result.i
        } else {
            result <- rbind(result, result.i)
        }
    }
    results[[j]] <- result
}

## creating graphs ------------------------------------------

main.cex=2.5
col1 <- "black"
pdf(file.path(FIG_DIR, filename),
    width=22,height=10) # for paper

par(mfrow=c(1,5))

for(p in 0:length(results)){

    if(p==0){
        par(cex.lab=1.3,
            mar=c(5,28,6,2.5),
            cex.lab=1.7, cex.axis=1.7,
            cex.main=main.cex)
        
        yloc <- c(23,22,21,18,17,16,13,12,11,8,7,6,3,2,1)
        pe <-  as.numeric(results[[1]][,1])
        plot(pe,yloc,
             ylim=c(1,24),
             xlim=c(0,0.001),
             xaxt='n',yaxt='n',bty='n',pch='',ylab='',xlab='')

        
        ## investment Provision
        mtext("Investment Protection", side=2, font=2,
              las=1, line=25, cex=1.75, at=24, adj=0, col=col1)
        mtext("Weak Legal Protection", side=2,
              las=1, line=15, cex=1.5, at=23, adj=0)
        mtext("Moderate Legal Protection", side=2,
              las=1, line=15, cex=1.5, at=22, adj=0)
        mtext("Strong Legal Protection", side=2,
              las=1, line=15, cex=1.5, at=21, adj=0)

        ## Reduction of tariff Barriers
        mtext("Reduction of Trade Barriers", side=2, font=2,
              las=1, line=25, cex=1.75, at=19, adj=0, col=col1)
        mtext("Low Reduction", side=2,
              las=1, line=15, cex=1.5, at=18, adj=0)
        mtext("Moderate Reduction", side=2,
              las=1, line=15, cex=1.5, at=17, adj=0)
        mtext("High Reduction", side=2,
              las=1, line=15, cex=1.5, at=16, adj=0)

        ## Export Subsidies
        mtext("Export Subsidies", side=2, font=2,
              las=1, line=25, cex=1.75, at=14, adj=0, col=col1)
        mtext("Low Subsidies", side=2,
              las=1, line=15, cex=1.5, at=13, adj=0)
        mtext("Moderate Subsidies", side=2,
              las=1, line=15, cex=1.5, at=12, adj=0)
        mtext("High Subsidies", side=2,
              las=1, line=15, cex=1.5, at=11, adj=0)

        ## Dispute Settlement
        mtext("Dispute Settlement Mechanism", side=2, font=2,
              las=1, line=25, cex=1.75, at=9, adj=0, col=col1)
        mtext("Passive Use", side=2,
              las=1, line=15, cex=1.5, at=8, adj=0)
        mtext("Moderate Use", side=2,
              las=1, line=15, cex=1.5, at=7, adj=0)
        mtext("Aggressive Use", side=2,
              las=1, line=15, cex=1.5, at=6, adj=0)

        ## Dispute Settlement
        mtext("Flexibility of Commitments", side=2, font=2,
              las=1, line=25, cex=1.75, at=4, adj=0, col=col1)
        mtext("Low Flexibility", side=2,
              las=1, line=15, cex=1.5, at=3, adj=0)
        mtext("Moderate Flxibiilty", side=2,
              las=1, line=15, cex=1.5, at=2, adj=0)
        mtext("High Flexibility", side=2,
              las=1, line=15, cex=1.5, at=1, adj=0)

    } else {
        par(cex.lab=1.3,
            mar=c(6,4,6,2.5),
            cex.lab=2, cex.axis=2,
            cex.main=main.cex)

        ## yloc <- seq(15, by=-1, length.out=15)
        yloc <- c(23,22,21,18,17,16,13,12,11,8,7,6,3,2,1)
        col1= "black"
        pe <-  as.numeric(results[[p]][,1])
        
        lo <- as.numeric(results[[p]][,1] - qnorm(.975)*results[[p]][,2])
        hi <- as.numeric(results[[p]][,1] + qnorm(.975)*results[[p]][,2])
        
        colors <- rep("black", length(pe))
        colors[c(3,12)] <- "red"

        plot(pe, yloc,
             xlim=c(-0.4, 0.5),
             ylim=c(1,24),
             pch=16,
             main=title[p],
             yaxt="n", ylab="", cex=2.5,
             col=colors,
             xlab="",
             )

        mtext("Change in Pr(Support Policy)",
              side=1, 
              line=4.5, cex=1.65)
        
        count <- 1
        for(i in 1:length(yloc)){
            loc <- yloc[i]
            segments(lo[count], loc, hi[count], loc, lwd=2, col=colors[i])
            count <- count + 1
        }
        abline(v=0, lty=2)
        xloc <- seq(-0.4,0.5,by=.1)
        for(k in xloc){
            abline(v=k, lty=2, col="grey80")
        }

        axis(side=2, at=yloc, tick=TRUE, labels=rep("",length(yloc)))
             

    }
}
dev.off()

